home *** CD-ROM | disk | FTP | other *** search
/ Libris Britannia 4 / science library(b).zip / science library(b) / DDJMAG / DDJ8607.ZIP / LETTERS.JUL < prev    next >
Text File  |  1986-07-31  |  8KB  |  284 lines
  1.  
  2.  
  3. ;Listing 1 - Square root algorithm plus code to benchmark and test it.
  4. ;
  5. INT_ROOT    SEGMENT
  6.         ASSUME  CS:INT_ROOT
  7. ;
  8. CALC_ROOT    PROC    NEAR
  9. ;----------------------------------------------;
  10. ; Argument passed in BX.                       ;
  11. ; Root returned in BX.                         ;
  12. ; All registers except BX preserved.           ;
  13. ;----------------------------------------------;
  14.         PUSH    AX
  15.         PUSH    CX
  16.         MOV    AX,BX        ;Hold argument in AX.
  17.         OR    BH,BH
  18.         JNZ    GE_256
  19.         CMP    BL,1        ;If arg is zero or one,
  20.         JBE    GOT_ROOT    ;then root = argument.
  21.         MOV    CL,4        ;If 2 <= argument <= 255
  22.         SHR    BL,CL        ;then guess = 3 + arg/16.
  23.         ADD    BL,3
  24.         JMP SHORT NEWTON
  25. GE_256:        MOV    BL,BH
  26.         MOV    BH,0
  27.         JS    GE_32768
  28.         CMP    BL,16
  29.         JAE    GE_4096
  30.         SHL    BL,1        ;If 256 <= argument <= 4095
  31.         SHL    BL,1        ;then guess=13+4*(arg hi byte).
  32.         ADD    BL,13
  33.         JMP SHORT NEWTON
  34. GE_4096:    ADD    BL,50        ;If 4096 <= argument <= 32767
  35.         JMP SHORT NEWTON    ;then guess = 50 + arg hi byte
  36. GE_32768:    CMP    BL,255     ;If arg hi byte=255 then root=255.
  37.         JZ    GOT_ROOT ;This prevents ovflow by DIV.
  38.         ADD    BL,40        ;If 32768 <= argument <= 65279
  39.         JNC    NEWTON        ;then guess = 40 + arg hi byte.    
  40.         MOV    BL,255        ;Guess must never exceed 255.
  41. ;
  42. NEWTON:        MOV    CX,AX        ;Save argument in CX.
  43.         DIV    BL        ;Divide by guess.
  44.         ADD    BL,AL        ;Guess + quotient.
  45.         RCR    BL,1        ;New guess=(old guess+quot)/2.
  46.         MOV    AL,BL         ;RCR shifts in carry from ADD.;
  47.         MUL    AL        ;If the square of the new guess
  48.         CMP    AX,CX        ;is greater than the argument,
  49.         JBE    GOT_ROOT    ;then we decrement new guess
  50.         DEC    BX        ;to get the correct root.
  51. ;
  52. GOT_ROOT:    POP    CX
  53.         POP    AX
  54.         RET
  55. ;
  56. CALC_ROOT    ENDP
  57. ;
  58. ;Listing 1 - Continued.
  59. ;------------------------------------------------------------;
  60. ; Code to time and test CALC_ROOT is designed to be run under;
  61. ; DEBUG and does not do a normal return to DOS but instead   ;
  62. ; does an INT 3 at the end of each routine.                  ;
  63. ;------------------------------------------------------------;
  64. TIME        PROC    FAR
  65. ;------------------------------------------------------------;
  66. ; TIME_ROOT computes the root of each of 65536 possible      ;
  67. ; arguments 15 times for a total of 983,040 roots.           ;
  68. ; TIME_OVER represents the looping overhead in TIME_ROOT.    ;
  69. ; The difference between the two times is the time to call   ;
  70. ; and execute CALC_ROOT.                                     ;
  71. ;------------------------------------------------------------;
  72.         MOV    BP,15
  73. TIME_OVER:    MOV    SI,0        ;Initial value.
  74. INNER_OVR:    MOV    BX,SI        
  75.         INC    SI
  76.         JNZ    INNER_OVR
  77.         DEC    BP
  78.         JNZ    TIME_OVER
  79. END_OVER:    MOV    AH,2
  80.         MOV    DL,7        ;Beep speaker.
  81.         INT    21H
  82.         INT     3
  83. ;
  84.         MOV    BP,15
  85. TIME_ROOT:    MOV    SI,0        ;Initial value.
  86. INNR_ROOT:    MOV    BX,SI
  87.         CALL    CALC_ROOT
  88.         INC    SI
  89.         JNZ    INNR_ROOT
  90.         DEC    BP
  91.         JNZ    TIME_ROOT
  92. END_TIME:    MOV    AH,2
  93.         MOV    DL,7        ;Beep speaker.
  94.         INT    21H
  95.         INT     3
  96. ;
  97. TIME        ENDP
  98. ;
  99. ;Listing 1 - Continued.
  100. ;
  101. TEST        PROC    FAR
  102. ;------------------------------------------------------------;
  103. ; If CHK_ROOT detects a bad root, it displays a message and  ;
  104. ; leaves the NG root in BX and the argument in SI.           ;
  105. ; If all roots are OK, a message to this effect is displayed.;
  106. ;------------------------------------------------------------;
  107.         MOV    SI,0        ;Initial value.
  108. CHK_ROOT:    MOV    BX,SI
  109.         CALL    CALC_ROOT
  110.         MOV    AX,BX
  111.         MUL    AX        ;DX,AX contains root^2.
  112.         JO    NG_ROOT        ;NG if root ^2 > 65535.
  113.         CMP    AX,SI
  114.         JA    NG_ROOT        ;NG if root^2 > argument.
  115.         MOV    AX,BX
  116.         INC    AX
  117.         MUL    AX        ;DX,AX contains (root+1)^2.
  118.         JO    NEXT_ARG    ;If ovflow then OK.
  119.         CMP    AX,SI
  120.         JBE    NG_ROOT        ;NG if (root+1)^2 <=argument.
  121. NEXT_ARG:    INC    SI
  122.         JNZ    CHK_ROOT
  123.         JMP SHORT OK_ROOTS
  124. ;
  125. OK_MSG    DB  0DH,0AH,'All roots tested OK.',0DH,0AH,'$'
  126. NG_MSG  DB  0DH,0AH,'Bad root in BX. Arg in SI.',0DH,0AH,'$'
  127. ;
  128. OK_ROOTS;    MOV    DX,OFFSET OK_MSG+100H      ;DS points to Pgm Seg
  129.         JMP SHORT DO_MSG        ;Pref which is 100H
  130. NG_ROOT:    MOV    DX,OFFSET NG_MSG+100H    ;lower than code seg.
  131. ;
  132. DO_MSG:        MOV    AH,9        ;Print string.
  133.         INT    21H
  134.         INT 3            ;Back to DEBUG.
  135. ;
  136. TEST        ENDP
  137. ;
  138. INT_ROOT    ENDS
  139. ;
  140.         END    TEST
  141.  
  142.  
  143.  
  144.  
  145. ;Listing 2 - BASIC program to test if a formula makes good square root guesses.
  146.  
  147. 10 FOR I = 2 TO 256
  148. 20 Q = I*I - 1
  149. 30 '
  150. 40 'Trial Formula to Calculate P0.
  151. 50 '
  152. 60 QHI = INT(Q/256): QLO = Q-QHI*256
  153. 70 IF QHI = 0 THEN P0 = INT(QLO/32) + 3: GOTO 160
  154. 80 IF QHI < 16 THEN P0 = 13 + 4*QHI: GOTO 160
  155. 90 IF QHI < 128 THEN P0 = QHI + 50; GOTO 160
  156. 100 IF QHI = 255 THEN P1 = 255: GOTO 210
  157. 110 P0 = QHI + 40
  158. 120 IF P0 > 255 THEN P0 = 255
  159. 130 '
  160. 140 ' Newtons Method
  161. 150 '
  162. 160 P1=INT((P0 + INT(Q / P0)) / 2)
  163. 170 IF P1 > 255 THEN PRINT "P1 > 255 when Q = ";Q:END
  164. 180 '
  165. 190 'Test result
  166. 200 '
  167. 210 P = INT(SQR(Q))
  168. 220 IF P1 <= P+1 GOTO 240
  169. 230 PRINT "For Q = ";Q;" P1 is greater than P+1. ":END
  170. 240 NEXT I
  171. 250 PRINT "Formula works for all worst cases."
  172.  
  173.  
  174.  
  175.  
  176. Listing Three
  177.  
  178.  INCLUDE MACLIB.ASM        ;by Neil R. Koozer
  179.  LIST ON            ;   Kellogg Star Rt. Box 125
  180.  MACLIST OFF            ;   Oakland, OR 97462
  181. ;                    (503)-459-3709
  182. ;SQR.ASM
  183.  
  184. ;Note that words like BR1 and BFS1 are macros to emulate BR:B and BFS:B
  185.  
  186.  GLOBAL SQR
  187.  
  188. SQR                ;square root function for 32000 floating point
  189.  RESTORE [R0]       ;use ret. addr as a pointer
  190.  MOVD 0(R0),R1     ;get operand address
  191.  MOVD 0(R1),R6     ;get part of operand
  192.  MOVD 4(R1),R7     ;get other part of operand
  193.  SBITB 31,R7       ;make the implicit 1 explicit
  194.  MOVD R6,R1        ;get exponent
  195.  ANDW 7FFH,R1      ;clean exponent
  196.  ADDW 3FFH,R1      ;fix offset
  197.  ROTW -1,R1        ;div exp by 2
  198.  CBITB 15,R1       ;test & clear wrap-around bit, 1=odd 0=even
  199.  SAVE [R0,R1]      ;save exp & ret addr
  200.  MOVB R7,R6        ;prepare for right shift
  201.  BFS1 SQR1         ;jump if exponent was odd
  202.  LSHD -4,R7        ;shift -3 for safety & -1 to get halfx
  203.  ROTD -4,R6
  204.  ANDW FF80H,R6     ;remove non-mantissa bits
  205.  MOVD 4C1BF828H,R3 ;y seed = 1.189.../2
  206.  BR1 SQR2
  207. SQR1
  208.  LSHD -3,R7        ;shift -3 for safety, -1 to get halfx, +1 because
  209.  ROTD -3,R6        ;  orig exp was odd
  210.  ANDW FF00H,R6     ;remove non-mantissa bits
  211.  MOVD 6A227E65H,R3 ;y seed = 1.68.../2
  212.  
  213.  
  214. SQR2               ;We will do 3 iterations with 32-bit precision
  215.  MOVD R7,R5        ;get halfx into R5
  216.  DEID R3,R4        ;R5 = halfx/y0  (the junk in R4 doesn't matter)
  217.  LSHD -1,R3        ;R3 = y0/2
  218.  ADDD R5,R3        ;R3 = new y0
  219.  
  220.  MOVD R7,R5        ;second iteration
  221.  DEID R3,R4
  222.  LSHD -1,R3
  223.  ADDD R5,R3
  224.  
  225.  MOVD R7,R5        ;third iteration
  226.  DEID R3,R4
  227.  LSHD -1,R3
  228.  ADDD R5,R3
  229.  
  230.  
  231.  MOVD R6,R4        ;Now the final iteration at full precision
  232.  MOVD R7,R5        ;get R4R5 = halfx from R6R7
  233.  
  234.  DEID R3,R4        ;now divide halfx (R4R5) by y (R2R3)
  235.  MOVD R2,R0
  236.  MEID R5,R0
  237.  NEGD R0,R0
  238.  SUBCD R1,R4
  239.  MOVD R4,R1
  240.  BCC1 DIV1
  241. DIV2
  242.  ADDQD -1,R5
  243.  ADDD R2,R0
  244.  ADDCD R3,R1
  245.  BCC DIV2
  246. DIV1
  247.  DEID R3,R0
  248.  MOVD R1,R4        ;R4R5 now = halfx / y
  249.  
  250.  MOVB R3,R2
  251.  LSHD -1,R3
  252.  ROTD -1,R2        ;R2R3 = y/2
  253.  
  254.  ADDD R4,R2
  255.  ADDCD R5,R3       ;R2R3 = y/2 + halfx/y
  256.                    ;4th iteration complete
  257.  
  258.  RESTORE [R0,R1]   ;get exponent & ret. addr
  259.  ADDD R2,R2        ;shift mantissa back where it belongs
  260.  ADDCD R3,R3
  261.  BCC1 SQR4         ;there should almost never be a carry
  262.  MOVB R3,R2
  263.  LSHD -1,R3        ;undo that last shift & zero the MSbit
  264.  ROTD -1,R2
  265.  ADDQW 1,R1        ;adjust exponent
  266.  BR1 SQR5          ;done
  267. SQR4
  268.  CBITB 31,R3       ;test & clear MSbit (make it a + sign)
  269.  BFS1 SQR5         ;it would virtually always be a 1
  270.  ADDD R2,R2        ;if not, shift left again
  271.  ADDCD R3,R3
  272.  ADDQW -1,R1       ;adjust exponent
  273.  CBITB 31,R3       ;make it +
  274. SQR5
  275.  ANDW F800H,R2     ;clean the mantissa
  276.  ORW R1,R2         ;append the exponent
  277.  MOVD 4(R0),R1     ;get addr of destination variable
  278.  MOVD R2,0(R1)     ;store result in dest. variable
  279.  MOVD R3,4(R1)
  280.  JUMP 8(R0)        ;return to caller
  281.  
  282.  END
  283.  
  284.